Protein and functional isoform levels and genetic variants of the BAFF and APRIL pathway components in systemic lupus erythematosus

Systemic lupus erythematosus (SLE) is the prototype of an autoimmune disease. Belimumab, a monoclonal antibody targets BAFF, is the only biologic approved for SLE and active lupus nephritis. BAFF is a cytokine with a key-regulatory role in the B cell homeostasis, which acts by binding to three receptors: BAFF-R, TACI and BCMA. TACI and BCMA also bind APRIL. Many studies reported elevated soluble BAFF and APRIL levels in the sera of SLE patients, but other questions about the role of this system in the disease remain open. The study aimed to investigate the utility of the cytokine levels in serum and urine as biomarkers, the role of non-functional isoforms, and the association of gene variants with the disease. This case–control study includes a cohort (women, 18–60 years old) of 100 patients (48% with nephritis) and 100 healthy controls. We used ELISA assays to measure the cytokine concentrations in serum (sBAFF and sAPRIL) and urine (uBAFF and uAPRIL); TaqMan Gene Expression Assays to quantify the relative mRNA expression of ΔBAFF, βAPRIL, and εAPRIL, and next-generation sequencing to genotype the cytokine (TNFSF13 and TNFSF13B) and receptor (TNFRSF13B, TNFRSF17 and TNFRSF13C) genes. The statistical tests used were: Kruskal–Wallis (qualitative variables), the Spearman Rho coefficient (correlations), the Chi-square and SKAT (association of common and rare genetic variants, respectively). As expected, sBAFF and sAPRIL levels were higher in patients than in controls (p ≤ 0.001) but found differences between patient subgroups. sBAFF and sAPRIL significantly correlated only in patients with nephritis (rs = 0.67, p ≤ 0.001) and βAPRIL levels were lower in patients with nephritis (p = 0.04), and ΔBAFF levels were lower in patients with dsDNA antibodies (p = 0.04). Rare variants of TNFSF13 and TNFRSF13B and TNFSF13 p.Gly67Arg and TNFRSF13B p.Val220Ala were associated with SLE. Our study supports differences among SLE patient subgroups with diverse clinical features in the BAFF/APRIL pathway. In addition, it suggests the involvement of genetic variants in the susceptibility to the disease.

Systemic lupus erythematosus (SLE) [OMIM #152700] is the prototype of autoimmune disease. It is a chronic disorder characterized by autoantibodies (auAb) production, immune complexes deposit, and heterogeneous clinical manifestations. Patients present a wide range of symptoms, including photosensitive rashes, discoid lesions, arthritis/arthralgia, nephritis, heart and lung alterations, and central nervous system disorders. Approximately 50-60% of patients with SLE have lupus nephropathy. This manifestation, especially in its most severe presentation, significantly increases disease-related morbidity and mortality 1,2 .
Belimumab is the only biologic approved for SLE and active lupus nephritis treatment by numerous regulatory agencies, including EMA and FDA. It is a monoclonal antibody that targets against B cell-activating factor belonging to the TNF family (BAFF, also called TNFSF13B, Blys, αTNF4, THANK and TALL-1), is a cytokine with a key-regulatory role in the B cell homeostasis. BAFF controls the number of peripheral B cells by binding to three receptors: BAFF-R (also called TNFRSF13C and BR3), TACI (or TNFRSF13B) and BCMA (also named BCM and TNFRSF17). BAFF-R only binds to BAFF and is essential for the survival and maturation of immature B cells. TACI, which has a higher affinity for a BAFF-like protein called APRIL (a proliferation-inducing ligand, also called TNFSF13), is a pivotal receptor in T cell-independent B response, in control of B cell-compartment size, and the isotype switching. BCMA has an intermediate affinity to BAFF or APRIL and promotes plasma cell survival. The signalling cascade through BAFF-R and BCMA stimulates the proliferation of B lymphocytes and counteracts apoptosis 3 . BAFF and APRIL are type II transmembrane proteins produced by myeloid linage cells, although lymphoid cells, including B and activated T cells, can also generate them 4 . BAFF becomes into the active soluble form after cleavage at the furin-protease site. APRIL mainly is produced in the soluble form previously by processing in the Golgi apparatus before releasing 5,6 . BAFF and APRIL, the same as the rest of the members of the TNF family, assemble as homotrimers. They can also join as heterotrimers, and the BAFF-APRIL recombinant-heterotrimers are biologically active 7,8 . The receptors are type III transmembrane proteins with similar gene and protein structures 3 . Many studies reported elevated levels of functional soluble BAFF and APRIL in the sera of SLE patients. Nevertheless, results regarding clinical features, activity and correlation between the levels of both cytokines are conflicting 9,10 . Very few studies have investigated the usefulness of urinary detection as a biomarker 11,12 . There are different isoforms of BAFF, products of alternative splicing, which could play a role in the pathogenesis of SLE. Especially interesting is ∆BAFF lacking exon 3, resulting in a shorter protein. Mouse ∆BAFF decreases bioactivity by associating with the complete forms become inactive heterotrimers. Similar events occur in APRIL, although with a non-well established functional impact. For instance, similarly to ∆BAFF, the βAPRIL isoform lacks an alternate in-frame exon in the central coding region, resulting in a shorter protein than canonical. In addition, no detection of sAPRIL in the supernatant of cells transfected with the β isoform has been reported 13 . εAPRIL isoform lacks exons 2 and 3 and represents a non-coding variant because the transcript is a candidate for nonsense-mediated mRNA decay.
SLE is a complex disease with a genetic basis and almost 30 genetic regions validated as predisposing to the disease 14 3 . Several studies reported the association of polymorphism in the genes of these cytokines and their receptors with autoimmune diseases, immunodeficiencies and cancer 15 .
To assess the involvement of the BAFF/APRIL system and its receptors in the development and clinical course of SLE, we performed a study evaluating the levels of the proteins in serum and urine and the role of non-functional isoforms, and the association of variants in cytokine and receptor genes.

Materials and methods
Study cohort. Quantification of the serum and urine levels of BAFF and APRIL. We used Enzyme-Linked Immunosorbent Assays (ELISA) to quantify the cytokine levels in serum and urine. In the case of BAFF (TNFSF13B), human in vitro Simple Step ELISA ® (ab188391-BAFF; Abcam, 330 Cambridge Science Park, Cambridge CB4 0FL, UK) was used according to the manufacturer's recommendations. Briefly, this system employed two antibodies, an affinity tag-labelled to capture and a reporter conjugated as the detector. The entire complex (capture antibody + analyte + detector antibody) is immobilized via immunoaffinity using an anti-tag antibody coated onto the well, testing 100 μL of each specimen (serum or urine) in each assay. The concentration, expressed in ng/mL, is calculated by optical density (OD) extrapolation with a standard curve (range 0-5 and the minimal detectable amount, MDA = 0.0127 ng/mL). APRIL (TNFSF13) Human in vitro ELISA kit (ab119505; Abcam) was employed according to the manufacturer's recommendations to quantify APRIL. This system uses APRIL specific antibodies on pre-coated plates to test 50 μL of each specimen (serum or urine) using an anti-APRIL biotinylated antibody as a second antibody. The OD is directly proportional to the APRIL amount captured on the plate. The concentration, expressed in ng/mL, is calculated by OD extrapolation with a standard curve (range 0-5 and MDA = 0.40 ng/mL).

Antinuclear antibodies screening and anti-dsDNA levels.
To investigate the presence of antinuclear antibodies (ANA), used ANA screen ELISA (Meridian Bioscience). To determine the specificities in positive sera, used the routinary method applied in the laboratory (EUROLINE ANA Profile 3 IgG kit, Euroimmun). To quantify anti-dsDNA antibody levels, tested a sample of 9 μL of serum from each participant with the Elia dsDNA test (Thermo Fisher Scientific) using a β-Galactosidase-conjugated secondary anti-IgG antibody. Anti-dsDNA antibody levels assignment, by comparing the fluorescence signals of the samples and calibrators (with known concentrations), and expressing the results in IU/mL (range 0.5-379 IU/mL) and classifying samples as neg < 10 IU/mL, doubtful 10-15 IU/mL and positive > 15 IU/mL. Quantification of mRNA levels of the BAFF and APRIL isoforms. For ΔBAFF and β and εAPRIL mRNA levels quantification, 10 7 peripheral blood mononuclear cells isolated obtained by density gradient were used to extract the total RNA with QIAmp RNA Mini Kits (Qiagen, Barcelona, Spain). We quantify total RNA by measurement of the OD260 and verify the integrity by electrophoresis and 260/280 nm absorption ratio. Concentrations ≥ 30 ng/μL and A260/280 absorbance ratios > 1.8 were considered acceptable, and stored samples were at − 80 °C until use.
The cDNA synthesis was performed by reverse-transcription using one µg of total RNA and random primers and the Superscript™ FirstStrand Synthesis System for RT-PCR (Thermo Fisher Scientific, Waltham, MA). www.nature.com/scientificreports/ Quantification of mRNA was performed by real-time PCR on a LightCycler 480 (Roche, Barcelona, Spain) using TaqMan Gene Expression Assays (ThermoFisher Scientific). (Supplementary Fig. 2). We tested the samples in duplicates, including a tube without any template (negative control) in each run. Data were analyzed with the LightCycler 4.05 software using the Calibrator Normalized Relative Quantification module with the efficiency correction method. A pool of cDNA from control samples was used as a calibrator to set a relative value of 1.
Only those values within the linear area of the standard curves were acceptable. Samples with Cp values > 35 and duplicates with a standard deviation of Cp > 0.3 were re-tested. The relative mRNA levels are expressed as the ratio of mRNA of the target isoform (ΔBAFF and β and εAPRIL) and the reference isoforms of the corresponding gene (BAFF isoform-1 and APRIL isoforms-α + γ) and normalized to the calibrator expression ratio.
Genotyping of the genes of the cytokines and their receptors. Following the manufacturer's instructions, we used peripheral blood samples to obtain DNA with QIamp DNA mini kits (Qiagen, Barcelona, Spain). Samples were stored at − 20 °C until used. DNA concentrations were quantified in a Qubit ® 3.0 fluorometer and diluted to 0.67 ng/µl. We genotype DNA samples by next-generation sequencing (NGS) in a custom-designed primer panel (AmpliSeq™ software, Ion Torrent, Thermo Fisher Scientific, Waltham, MA). This panel targeted all the coding regions and flanking intronic sequences of the five genes included in the study: . Quantitative values display as median, interquartile range (IQR, Q3-Q1). The variable-correlation study includes only individuals without missing data, using the Spearman Rho coefficient because data exhibits a non-linear association between variables and contains outliers (https:// www. socsc istat istics. com/ tests/ spear man/). When significant, the r s -values 0-0.20 were considered very weak correlation, 0.20-0.39 weak correlation, 0.40-0.69 moderate correlation, 0.70-0.89 strong correlation, and 0.90-1.00 very strong correlation. In all the comparisons, the p-values < 0.05 were considered significant.
In the association study, minor allele frequencies (MAF) of genetic variants found in our study (91 patients and 91 controls) were obtained through the Cellbase database 19 using the 1000 Genomes Project data source. This way, variants were classified as commons (MAF ≥ 0.05), rares (MAF < 0.05) or undetected (MAF = 0). We catalogued all the variants according to the American College of Medical Genetics and Genomics (ACMG) criteria using Varsome (https:// varso me. com/) as Benign (B), Likely Benign (LB), Variant of Uncertain Significance (VUS), likely pathogenic (LP), and pathogenic (P). Common and rare variants were treated separately in the case-control study. We used the Chi-square to evaluate the association of the common variants, and p-values were adjusted using the false discovery rate (FDRadj.pval). To analyze the association of the rare variants, we used Vartools, a module of the Zhbannikov software (http:// izhba nnikov. github. io/ varto ols/) and applied three different strategies. The first included all the rare variants in each gene (MAF ≤ 0.05 in 1 KG-ALL). The second, those rare variants with possible functional relevance (MAF ≤ 0.05 in 1 KG-ALL annotated as missense and frameshift variants). Lastly, the third included those rare variants with a probable clinical significance (MAF ≤ 0.05 catalogued as VUS, LP or P according to the ACMG criteria). The analysis with the SNP-set Kernel Association Test (SKAT) 20 includes those genes with at least two rare variants that meet the filter criteria. A p-value < 0.05 was considered significant.
Ethical approval. The Hospitales Universitarios Virgen Macarena and Virgen del Rocio's ethical committee approved this study (peiba_DictamenFavorable2019114215856). The authors declared that all methods were carried out following relevant guidelines and regulations. Table 1 displays the characteristics of the cohort. Patients were selected so that the proportion of those with and without nephritis was approximately 1:1. Regarding the autoantibodies, 85% of patients had a positive result in the ANA test in the serum obtained for this study, and only 27% had anti-dsDNA antibodies in the same serum. The most common autoantibodies specificities were Ro/La (27%) and Sm/RNP (19%). Most of the patients had a non-active disease and received treatment with hydroxychloroquine in monotherapy or combination.

Quantification of the serum and urine levels of BAFF and APRIL. Median levels of BAFF and
APRIL in serum (sBAFF and sAPRIL) were higher in patients than in controls (p < 0.001 in both cases) (Fig. 1). Both serum cytokine levels significantly correlated in patients with nephritis but not in patients without or in controls ( Table 2). There were no significant differences in the median levels of both cytokines between patients with and without nephritis ( Supplementary Fig. 3 Only two patients (with dsDNA antibodies, without nephritis and with a mild activity) had sBAFF concentration values over 2.0 ng/mL. There were no significant differences in the features of the groups of patients with sAPRIL values greater or less than 4.0 (N = 15: 3, 20%, with anti-dsDNA antibodies, 7, 47% with nephritis, 5, 33% with active disease vs. N = 72: 21, 29% with anti-dsDNA antibodies, 35, 49% with nephritis, 26, 36% with active disease).

Quantification of mRNA levels of the ΔBAFF and β and εAPRIL isoforms.
Concerning the relative expression of non-functional isoforms in patients and controls, there were no significant differences in the cases of ΔBAFF and βAPRIL However, the εAPRIL median relative expression was lower in patients than in controls (Fig. 2a). For patient subgroups, the relative expression of βAPRIL was lower in patients with nephritis than in patients without (Fig. 2b), and the relative levels of ΔBAFF were lower in patients with dsDNA antibodies than in patients without (Fig. 2c). The rest of the comparisons were not statistically significant (data not shown).
The relative expression levels of β and εAPRIL correlated in patients and controls, whereas the relative expression levels of εAPRIL with levels of sAPRIL only correlated in controls. Lastly, the relative expression levels of the βAPRIL and ΔBAFF did not correlate with the serum levels of the corresponding protein in any group (Table 2). Table 1 displayed detailed information on the 52 variants found in the five genes included in the study. Table 3 Fig. 6).

Association study of the genes of the cytokines and their receptors. Supplementary
From a total of 16 common variants (MAF > 0.05 in the 1000 Genome all population Database), only in one case, TNFSF13 p.Gly67Arg (rs11552708), the distribution was significantly different in patients and controls (MAF in patients 0.10, MAF in controls 0.03 p = 0.014, OR 3.21 95% CI 1.25-9.02), although the FDRadj. pval > 0.05 (Table 4).  www.nature.com/scientificreports/ Regarding the rare variants, TNFRSF13B Val220Ala (considered B) was the only associated according to an individual test (MAF in patients 0.0 vs. 0.02 in controls, p = 0.04. FDRadj.pval > 0.05). Concerning the analysis of rare variants at the gene level, TNFRSF13B was associated with SLE using SKAT with all the three filters applied: including all the variants with MAF < 0.05 (p = 0.03), restricting the analysis to those rare variants missense or frameshift (p = 0.006) and including only rare variants VUS, LP and P (p = 0.01). The TNFRSF13B rare variant number was lower in patients (20/182 alleles) than in controls (37/182 alleles). Also, the number of individuals with TNFRSF13B rare variants was lesser in patients (15/91, 17%) than in controls (24/91, 26%), although the difference was not statistically significant (p = 0.052). TNFSF13B was significantly associated with SKAT with rare variants VUS, LP and P (p = 0.04) though it did not reach significance with all the rare variants included in the test (p = 0.08). For the remaining genes, statistical significance was lacking. TNFSF13 was near statistical  www.nature.com/scientificreports/ significance with SKAT in the filter with all the rare variants included in the test (p = 0.057), but this gene did not meet the inclusion criteria with the other two filters (Table 5).

Discussion
In agreement with most of those previously published, our study reports higher levels of sBAFF and sAPRIL in female SLE patients than in healthy controls ethnically, gender and age-matched (revised in 9,10 ). Our results suggest no influence of the serum levels of these cytokines with the dsDNA antibodies or the activity. Several studies reported discrepancies regarding the relationship between the concentration of these cytokines, the clinical features and the activity score 9,10 . Characteristics of the different cohorts and the index of activity used may influence the results. In this sense, our cohort consisted of patients with a low activity index. Also, the number of patients with and without nephritis is equivalent, but it is unbalanced regarding the dsDNA antibodies and the activity. These characteristics could influence the relatively low serum concentration of both cytokines found.
In the case of sBAFF, only two patients had concentrations over those established as an independent prognostic  Table 5. Analysis of the association of rare variants in the five genes included in the study. 1 SNP-set (Sequence) Kernel Association Test Method (SKAT, Asymptotic p-value). This test analyzes the association of rare variants (Minor Allele Frequency, MAF < 0.05) grouped by the gene. The p-values < 0.05 were considered significant (shown in bold). VUS variant of uncertain significance, LP likely pathogenic, P pathogenic. Only the genes that met the evaluation criteria are displayed. The rest of the genes could not be evaluated with the corresponding filter because they did not have at least two rare variants that meet the condition. www.nature.com/scientificreports/ factor for flares (2 ng/mL) 21 . Regarding sAPRIL, patients with sAPRIL values above and below the cut-off point reported as a predictor of the resistance to treatment (4.0 ng/mL) 22 were not significantly different concerning the clinical variables analyzed. A previously investigated question that has reported contradictory results is the coexistence of both molecules in the serum of SLE patients [23][24][25] . Our results strongly support the coexistence of both cytokines in the serum of patients with nephritis with a moderate correlation. This correlation in the subgroup of patients with nephritis may be due to parallel changes in the levels of both cytokines, but also to the fact that the levels of heterotrimers are higher in this subgroup (because both BAFF and APRIL ELISA assays detect heterotrimers) Previous studies reported higher levels of heterotrimers in patients with SLE compared with healthy control and patients with rheumatoid arthritis 7,8 . The restriction to patients with nephritis could explain the contradictory results on the coexistence of both molecules in SLE patient sera and be relevant in the treatment options. Thus, Belimumab, which blocks BAFF without neutralizing APRIL, has a minimal inhibitory activity on heterotrimers BAFF 2 APRIL 1 and no activity in BAFF 1 APRIL 2 26,27 . Theoretically, patients with the most elevated levels of both cytokines or heterotrimers would respond worse to the Belimumab treatment. There are studies which are reported differences in the kinetic sBAFF and sAPRIL in immunotherapy without BAFF blockade. The sAPRIL levels decreased after treatment predicting the response in patients with proliferative lupus nephritis. sBAFF levels under 1.5 ng/mL at baseline are a good predictor of response in patients with lupus nephritis but remain unchanged afterwards. So, even with conventional immunotherapy, sBAFF and sAPRIL could be used as biomarkers predictors of response 28 . Further studies are necessary to establish how the coexistence of both cytokines or levels of heterotrimers influences the response to treatments.
Another remaining question is about urinary excretion of BAFF and APRIL, proposed as biomarkers of SLE. The urinary excretion of BAFF seems to be an uncommon event, whereas the urinary excretion of APRIL may be a physiological process because it was detected often in controls. Although these results may seem surprising, similar findings have previously been published 29 . Therefore, our results do not support the utility of these urinary cytokines as biomarkers of SLE. It hypothesized that the urinary excretion of the cytokines could reduce their serum levels. In our study, the correlation between uAPRIL and sAPRIL levels in controls is robust, although moderate, and lacking in patients. The meaning of the lack of correlation between uAPRIL and sAPRIL levels in patients could be related to the disease though this result would need replication.
BAFF and APRIL, as a result of alternative splicing events, have different mature forms of mRNA. The nonfunctional isoforms such as ΔBAFF, βAPRIL and εAPRIL may regulate the active form of the protein, but their functional impact is not well characterized yet. To our knowledge, there are no previous assessments of the role of the non-functional isoforms of these cytokines in SLE. The mRNA simultaneous transcription of the two isoforms of APRIL may be a physiological process because their levels correlate in patients and controls. εAPRIL, underproduced in the patient group, is a non-coding variant which could have regulatory functions 30 . The positive correlation between εAPRIL and sAPRIL in healthy controls that is lacking in patients could indicate that, physiologically, the production of this isoform is a mechanism of response to elevated protein levels. This control mechanism may be disturbed in patients, the group in which this isoform is underproduced. Our results do not support any influence of βAPRIL and ΔBAFF on the serum protein levels because there was no correlation in patient or control groups. In addition, βAPRIL was underproduced in patients with nephritis and ΔBAFF in patients with anti-dsDNA antibodies, suggesting an influence of the regulatory activity of these isoforms in some patient subgroups.
Many studies investigated the genetic factors related to SLE by comparing the distribution of single nucleotide polymorphisms (SNPs) in patients and controls. The high-throughput genetic studies (such as Genome-Wide Association Studies, GWAS) have been very informative in SLE and other polygenic diseases. Nevertheless, GWAS are not approaches designed to detect rare variants, which are lost even with imputation processes because their extremely-low frequency. The NGS permits the analysis of the common and the rare variants. This finemapping approach could help solve the missing heritability problem and validate therapeutic targets. Concerning the common variants included in our study, the results suggest an association of the TNFSF13 p.Gly67Arg in our population (risk variant Arg67), although the p-value became non-significant after correction. Several studies reported an association of this position with SLE in Japanese and, to a lesser extent, in African-Americans and Hispanics (risk variant Gly67), but not in American white people [31][32][33] . TNFSF13 p.Gly67Arg is a B variant according to the ACMG criteria, and its functional significance remains unknown. Altogether, results suggest an association with the region but rule out p.Gly67Arg as the causal variant. In agreement with most of the previously published studies [34][35][36][37] , there were no common variants associated with SLE in the other four genes inside the region included in our study [34][35][36][37][38] , although the association with other diseases has been described [39][40][41][42] . Regarding the individual analysis of rare variants, TNFRSF13B p.Val220Ala was associated with the disease, being Ala220 lesser common in patients. This variant is considered B, but it could disturb the stability centres of the protein 42 . In any case, the individual associations study of rare variants in complex diseases has difficulties because of the lack of statistical power. To get around this question, specific statistical methods: grouping variants and checking the association of the entire gene, have been designed. In the present study, TNFRSF13B was associated with the disease independently of the strategy used. Variants in TNFRSF13B have been associated with antibody deficiency, finding a high frequency of variants in patients with common variable immunodeficiency and IgA deficiency 43 . But, the accumulation of variants in TNFRSF13B could confer some evolutionary advantage since the gene has an unexpected diversity and the IgA deficiency is often asymptomatic 44 . In this sense, the number of rare variants in TNFRSF13B in our cohort of patients was lower than in controls. The other gene with associated rare variants was TNFSF13B, but only in the analysis restricted to variants VUS + LP + P. Previous studies reported a rare INDEL variant GCTGT>A located outside the region included in our study (in the 3´UTR region of TNFSF13B gene) as associated with SLE and other autoimmune diseases 45  www.nature.com/scientificreports/ the remaining genes, although their study has more limitations and, in general, our findings with rare variants need replication in other cohorts.
In conclusion, our study supports differences among patient subgroups in the coexistence of both cytokines or the levels of heterotrimers. It suggests a role of the non-functional isoforms that may also be related to clinical features. In addition, it supports the involvement of rare variants in TNFSF13B and TNFRSF13B in the disease.

Data availability
All data relevant to the study are included in the article or uploaded as supplementary information. This study is registered in ClinicalTrials.gov (https:// regis ter. clini caltr ials. gov/ prs) NCT 03919643 (Initial release 02/14/2019). Further data can be made available upon request mariaf.gonzalez.sspa@juntadeandalucia.es.